#Set-up
rm(list=ls())
library(coda)
library(lattice)
library(rjags)
library(xtable)
library(gmodels)
#setwd('/Users/monogan/Documents/looseData/coda2016/')
#setwd('/Users/jamie/looseData/coda2016/')
#setwd('C:/Users/monogan/Desktop/coda2016')
#setwd("/Volumes/MONOGAN/psrsd/anxiety2016/data/REPLICATION/")
#setwd("/Users/monogan/Documents/psrsd/anxiety2016/data/REPLICATION/")


### 2016 MODEL ###
data<-read.table("anes16.txt",header=T)

#descriptives
data.desc<-subset(data,select=-caseid)
print(cbind(
apply(data.desc,2,mean),
apply(data.desc,2,sd),
apply(data.desc,2,min),
apply(data.desc,2,max)
),digits=3)

Sys.time()
#define model with default burn-in of 1000 iterations 
mod<-jags.model(file="MODELonly2016.txt", data=list(N=dim(data)[1],
	vote.rep=data$vote.rep,pid_x=data$pid_x,issues=data$issues,personal=data$personal,
	ang.own=data$ang.own,hp.own=data$hp.own,afr.own=data$afr.own,prd.own=data$prd.own,disg.own=data$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod,n.iter=49000) 
#update(mod,n.iter=1000) 

#post burn-in sampling
output.mod<-coda.samples(mod,c("beta","constant","alpha"),500000)
#output.mod<-coda.samples(mod,c("beta","constant","alpha"),1000)

#Density Plots
pdf("anxPid2016.pdf",pointsize=8,width=4,height=4)
anx.pid<-c(output.mod[[1]][,"beta[5]"],output.mod[[2]][,"beta[5]"],output.mod[[3]][,"beta[5]"])
plot(density(anx.pid),xlab="Coefficient",main="Interaction Between Anxiety and Partisanship")
dev.off()

pdf("anxIss2016.pdf",pointsize=8,width=4,height=4)
anx.iss<-c(output.mod[[1]][,"beta[6]"],output.mod[[2]][,"beta[6]"],output.mod[[3]][,"beta[6]"])
plot(density(anx.iss),xlab="Coefficient",main="Interaction Between Anxiety and Issue Distance")
dev.off()

pdf("anxPers2016.pdf",pointsize=8,width=4,height=4)
anx.pers<-c(output.mod[[1]][,"beta[7]"],output.mod[[2]][,"beta[7]"],output.mod[[3]][,"beta[7]"])
plot(density(anx.pers),xlab="Coefficient",main="Interaction Between Anxiety and Candidate Qualities")
dev.off()

#compute DIC
dic.samples(mod,10000)

#obtain model output
summary.cumulative<-summary(output.mod,quantiles=c(0.05,0.5,0.95))
summary.cumulative$statistics
summary.cumulative$quantiles
xtable(cbind(summary.cumulative$statistics[,1:2],summary.cumulative$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod)
gelman.diag(output.mod)
Sys.time()

###Predicted Probabilities 2016###
#Anxiety is on a logit 0 to 1 scale
plogis(0)
plogis(1)
plogis(-1)

#bind samples into a matrix
output.mat<-rbind(output.mod[[1]],output.mod[[2]],output.mod[[3]])
output.mat<-output.mat[,-c(1:5)]
output.means<-apply(output.mat,2,mean)

#create vectors of average predictor values at low, medium, and high anxiety
low.averages<-c(0.482,0.509,0.518,0.2689414,0.2689414*0.482,0.2689414*0.509,0.2689414*0.518,1)
med.averages<-c(0.482,0.509,0.518,0.5,0.5*0.482,0.5*0.509,0.5*0.518,1)
high.averages<-c(0.482,0.509,0.518,0.7310586,0.7310586*0.482,0.7310586*0.509,0.7310586*0.518,1)
###plogis(low.averages%*%output.means)
###plogis(med.averages%*%output.means)
###plogis(high.averages%*%output.means)

##create matrices where PARTY varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.party<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.party[,1]<-ruler
low.averages.party[,5]<-ruler*0.2689414
med.averages.party<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.party[,1]<-ruler
med.averages.party[,5]<-ruler*0.5
high.averages.party<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.party[,1]<-ruler
high.averages.party[,5]<-ruler*0.7310586

#get averages
low.averages.party.prob<-plogis(low.averages.party%*%output.means)
med.averages.party.prob<-plogis(med.averages.party%*%output.means)
high.averages.party.prob<-plogis(high.averages.party%*%output.means)

#forecast for each iteration
low.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.party.all[i,]<-plogis(low.averages.party%*%output.mat[i,])
	}
low.party.ci.low<-apply(low.party.all,2,quantile,.05)
low.party.ci.high<-apply(low.party.all,2,quantile,.95)

med.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.party.all[i,]<-plogis(med.averages.party%*%output.mat[i,])
	}
med.party.ci.low<-apply(med.party.all,2,quantile,.05)
med.party.ci.high<-apply(med.party.all,2,quantile,.95)

high.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.party.all[i,]<-plogis(high.averages.party%*%output.mat[i,])
	}
high.party.ci.low<-apply(high.party.all,2,quantile,.05)
high.party.ci.high<-apply(high.party.all,2,quantile,.95)

pdf("partyProb2016.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.party.ci.low,rev(med.party.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.party.ci.low,rev(low.party.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.party.ci.low,rev(high.party.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("partyProb2016mean.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where ISSUE proximity varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.issue<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.issue[,2]<-ruler
low.averages.issue[,6]<-ruler*0.2689414
med.averages.issue<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.issue[,2]<-ruler
med.averages.issue[,6]<-ruler*0.5
high.averages.issue<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.issue[,2]<-ruler
high.averages.issue[,6]<-ruler*0.7310586

#get averages
low.averages.issue.prob<-plogis(low.averages.issue%*%output.means)
med.averages.issue.prob<-plogis(med.averages.issue%*%output.means)
high.averages.issue.prob<-plogis(high.averages.issue%*%output.means)

#forecast for each iteration
low.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.issue.all[i,]<-plogis(low.averages.issue%*%output.mat[i,])
	}
low.issue.ci.low<-apply(low.issue.all,2,quantile,.05)
low.issue.ci.high<-apply(low.issue.all,2,quantile,.95)

med.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.issue.all[i,]<-plogis(med.averages.issue%*%output.mat[i,])
	}
med.issue.ci.low<-apply(med.issue.all,2,quantile,.05)
med.issue.ci.high<-apply(med.issue.all,2,quantile,.95)

high.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.issue.all[i,]<-plogis(high.averages.issue%*%output.mat[i,])
	}
high.issue.ci.low<-apply(high.issue.all,2,quantile,.05)
high.issue.ci.high<-apply(high.issue.all,2,quantile,.95)

pdf("issueProb2016.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Trump Issue Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.issue.ci.low,rev(med.issue.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.issue.ci.low,rev(low.issue.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.issue.ci.low,rev(high.issue.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("issueProb2016mean.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Trump Issue Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where PERSONAL qualities vary for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.personal<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.personal[,3]<-ruler
low.averages.personal[,7]<-ruler*0.2689414
med.averages.personal<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.personal[,3]<-ruler
med.averages.personal[,7]<-ruler*0.5
high.averages.personal<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.personal[,3]<-ruler
high.averages.personal[,7]<-ruler*0.7310586

#get averages
low.averages.personal.prob<-plogis(low.averages.personal%*%output.means)
med.averages.personal.prob<-plogis(med.averages.personal%*%output.means)
high.averages.personal.prob<-plogis(high.averages.personal%*%output.means)

#forecast for each iteration
low.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.personal.all[i,]<-plogis(low.averages.personal%*%output.mat[i,])
	}
low.personal.ci.low<-apply(low.personal.all,2,quantile,.05)
low.personal.ci.high<-apply(low.personal.all,2,quantile,.95)

med.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.personal.all[i,]<-plogis(med.averages.personal%*%output.mat[i,])
	}
med.personal.ci.low<-apply(med.personal.all,2,quantile,.05)
med.personal.ci.high<-apply(med.personal.all,2,quantile,.95)

high.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.personal.all[i,]<-plogis(high.averages.personal%*%output.mat[i,])
	}
high.personal.ci.low<-apply(high.personal.all,2,quantile,.05)
high.personal.ci.high<-apply(high.personal.all,2,quantile,.95)

pdf("personalProb2016.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Trump Personal Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.personal.ci.low,rev(med.personal.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.personal.ci.low,rev(low.personal.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.personal.ci.low,rev(high.personal.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("personalProb2016mean.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Trump Personal Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

###Alternate Predicted Probabilities 2016, setting anxiety at empirical minimum, mean, and maximum.###
#Anxiety is on a logit 0 to 1 scale
plogis(0)#now .5018
plogis(1)#now .9162
plogis(-1)# now .1981

#bind samples into a matrix
output.mat<-rbind(output.mod[[1]],output.mod[[2]],output.mod[[3]])
output.mat<-output.mat[,-c(1:5)]
output.means<-apply(output.mat,2,mean)

#create vectors of average predictor values at low, medium, and high anxiety
low.averages<-c(0.482,0.509,0.518,.1981,.1981*0.482,.1981*0.509,.1981*0.518,1)
med.averages<-c(0.482,0.509,0.518,.5018,.5018*0.482,.5018*0.509,.5018*0.518,1)
high.averages<-c(0.482,0.509,0.518,.9162,.9162*0.482,.9162*0.509,.9162*0.518,1)
###plogis(low.averages%*%output.means)
###plogis(med.averages%*%output.means)
###plogis(high.averages%*%output.means)

##create matrices where PARTY varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.party<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.party[,1]<-ruler
low.averages.party[,5]<-ruler*.1981
med.averages.party<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.party[,1]<-ruler
med.averages.party[,5]<-ruler*0.5018
high.averages.party<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.party[,1]<-ruler
high.averages.party[,5]<-ruler*.9162

#get averages
low.averages.party.prob<-plogis(low.averages.party%*%output.means)
med.averages.party.prob<-plogis(med.averages.party%*%output.means)
high.averages.party.prob<-plogis(high.averages.party%*%output.means)

#forecast for each iteration
low.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.party.all[i,]<-plogis(low.averages.party%*%output.mat[i,])
	}
low.party.ci.low<-apply(low.party.all,2,quantile,.05)
low.party.ci.high<-apply(low.party.all,2,quantile,.95)

med.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.party.all[i,]<-plogis(med.averages.party%*%output.mat[i,])
	}
med.party.ci.low<-apply(med.party.all,2,quantile,.05)
med.party.ci.high<-apply(med.party.all,2,quantile,.95)

high.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.party.all[i,]<-plogis(high.averages.party%*%output.mat[i,])
	}
high.party.ci.low<-apply(high.party.all,2,quantile,.05)
high.party.ci.high<-apply(high.party.all,2,quantile,.95)

pdf("partyProb2016Alt.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.party.ci.low,rev(med.party.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.party.ci.low,rev(low.party.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.party.ci.low,rev(high.party.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("partyProb2016meanAlt.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where ISSUE proximity varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.issue<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.issue[,2]<-ruler
low.averages.issue[,6]<-ruler*.1981
med.averages.issue<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.issue[,2]<-ruler
med.averages.issue[,6]<-ruler*0.5018
high.averages.issue<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.issue[,2]<-ruler
high.averages.issue[,6]<-ruler*.9162

#get averages
low.averages.issue.prob<-plogis(low.averages.issue%*%output.means)
med.averages.issue.prob<-plogis(med.averages.issue%*%output.means)
high.averages.issue.prob<-plogis(high.averages.issue%*%output.means)

#forecast for each iteration
low.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.issue.all[i,]<-plogis(low.averages.issue%*%output.mat[i,])
	}
low.issue.ci.low<-apply(low.issue.all,2,quantile,.05)
low.issue.ci.high<-apply(low.issue.all,2,quantile,.95)

med.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.issue.all[i,]<-plogis(med.averages.issue%*%output.mat[i,])
	}
med.issue.ci.low<-apply(med.issue.all,2,quantile,.05)
med.issue.ci.high<-apply(med.issue.all,2,quantile,.95)

high.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.issue.all[i,]<-plogis(high.averages.issue%*%output.mat[i,])
	}
high.issue.ci.low<-apply(high.issue.all,2,quantile,.05)
high.issue.ci.high<-apply(high.issue.all,2,quantile,.95)

pdf("issueProb2016Alt.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Trump Issue Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.issue.ci.low,rev(med.issue.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.issue.ci.low,rev(low.issue.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.issue.ci.low,rev(high.issue.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("issueProb2016meanAlt.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Trump Issue Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where PERSONAL qualities vary for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.personal<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.personal[,3]<-ruler
low.averages.personal[,7]<-ruler*.1981
med.averages.personal<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.personal[,3]<-ruler
med.averages.personal[,7]<-ruler*0.5018
high.averages.personal<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.personal[,3]<-ruler
high.averages.personal[,7]<-ruler*.9162

#get averages
low.averages.personal.prob<-plogis(low.averages.personal%*%output.means)
med.averages.personal.prob<-plogis(med.averages.personal%*%output.means)
high.averages.personal.prob<-plogis(high.averages.personal%*%output.means)

#forecast for each iteration
low.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.personal.all[i,]<-plogis(low.averages.personal%*%output.mat[i,])
	}
low.personal.ci.low<-apply(low.personal.all,2,quantile,.05)
low.personal.ci.high<-apply(low.personal.all,2,quantile,.95)

med.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.personal.all[i,]<-plogis(med.averages.personal%*%output.mat[i,])
	}
med.personal.ci.low<-apply(med.personal.all,2,quantile,.05)
med.personal.ci.high<-apply(med.personal.all,2,quantile,.95)

high.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.personal.all[i,]<-plogis(high.averages.personal%*%output.mat[i,])
	}
high.personal.ci.low<-apply(high.personal.all,2,quantile,.05)
high.personal.ci.high<-apply(high.personal.all,2,quantile,.95)

pdf("personalProb2016Alt.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Trump Personal Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.personal.ci.low,rev(med.personal.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.personal.ci.low,rev(low.personal.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.personal.ci.low,rev(high.personal.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("personalProb2016meanAlt.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Trump Personal Advantage",ylab="Probability of Trump Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()



### MODEL WITHOUT DISGUST IN 2016 ###
data<-read.table("anes16.txt",header=T)

Sys.time()
#define model with default burn-in of 1000 iterations 
mod<-jags.model(file="MODELannual.txt", data=list(N=dim(data)[1],
	vote.rep=data$vote.rep,pid_x=data$pid_x,issues=data$issues,personal=data$personal,
	ang.own=data$ang.own,hp.own=data$hp.own,afr.own=data$afr.own,prd.own=data$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod,n.iter=49000) 

#post burn-in sampling
output.mod<-coda.samples(mod,c("beta","constant","alpha"),500000)

#compute DIC
dic.samples(mod,10000)

#obtain model output
summary.cumulative<-summary(output.mod,quantiles=c(0.05,0.5,0.95))
summary.cumulative$statistics
summary.cumulative$quantiles
xtable(cbind(summary.cumulative$statistics[,1:2],summary.cumulative$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod)
gelman.diag(output.mod)
Sys.time()



### SELECTION MODEL IN 2016 ###
rm(list=ls())
selection<-read.table("anes16selection.txt",header=T)
#selection<-read.table("anes16.txt",header=T)

#define model with default burn-in of 1000 iterations 
mod.ideol<-jags.model(file="MODELheckman2016.txt", data=list(N=dim(selection)[1],start=341,
	vote.rep=selection$vote.rep,pid_x=selection$pid_x,issues=selection$issues,personal=selection$personal,
	ang.own=selection$ang.own,hp.own=selection$hp.own,afr.own=selection$afr.own,prd.own=selection$prd.own,
	disg.own=selection$disg.own,voted=selection$voted),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.ideol,n.iter=49000) 
#update(mod.ideol,n.iter=1000) 

#post burn-in sampling
output.mod.ideol<-coda.samples(mod.ideol,c("beta","constant","alpha","gamma","intercept","rho"),100000)
#output.mod.ideol<-coda.samples(mod.ideol,c("beta","constant","alpha","gamma","intercept","rho"),1000)

#compute DIC
dic.samples(mod.ideol,10000)

#obtain model output
summary.cumulative.ideol<-summary(output.mod.ideol,quantiles=c(0.05,0.5,0.95))
summary.cumulative.ideol$statistics
summary.cumulative.ideol$quantiles
xtable(cbind(summary.cumulative.ideol$statistics[,1:2],summary.cumulative.ideol$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.ideol)
gelman.diag(output.mod.ideol)
Sys.time()




### THERMOMETER IN 2016 ###
rm(list=ls())
thermometer<-read.table("thermData.txt",header=T)

therm.cor<-matrix(NA,nrow=10,ncol=2)
therm.cor[1,1]<-cor(thermometer$dAngry,thermometer$demTherm)
therm.cor[2,1]<-cor(thermometer$dAfraid,thermometer$demTherm)
therm.cor[3,1]<-cor(thermometer$dDisgust,thermometer$demTherm)
therm.cor[4,1]<-cor(thermometer$dHope,thermometer$demTherm)
therm.cor[5,1]<-cor(thermometer$dProud,thermometer$demTherm)

therm.cor[6,1]<-cor(thermometer$rAngry,thermometer$demTherm)
therm.cor[7,1]<-cor(thermometer$rAfraid,thermometer$demTherm)
therm.cor[8,1]<-cor(thermometer$rDisgust,thermometer$demTherm)
therm.cor[9,1]<-cor(thermometer$rHope,thermometer$demTherm)
therm.cor[10,1]<-cor(thermometer$rProud,thermometer$demTherm)

therm.cor[1,2]<-cor(thermometer$dAngry,thermometer$repTherm)
therm.cor[2,2]<-cor(thermometer$dAfraid,thermometer$repTherm)
therm.cor[3,2]<-cor(thermometer$dDisgust,thermometer$repTherm)
therm.cor[4,2]<-cor(thermometer$dHope,thermometer$repTherm)
therm.cor[5,2]<-cor(thermometer$dProud,thermometer$repTherm)

therm.cor[6,2]<-cor(thermometer$rAngry,thermometer$repTherm)
therm.cor[7,2]<-cor(thermometer$rAfraid,thermometer$repTherm)
therm.cor[8,2]<-cor(thermometer$rDisgust,thermometer$repTherm)
therm.cor[9,2]<-cor(thermometer$rHope,thermometer$repTherm)
therm.cor[10,2]<-cor(thermometer$rProud,thermometer$repTherm)

colnames(therm.cor)<-c("Clinton Thermometer","Trump Thermometer")
rownames(therm.cor)<-c("Clinton Angry","Clinton Afraid","Clinton Disgust","Clinton Hope","Clinton Proud","Trump Angry","Trump Afraid","Trump Disgust","Trump Hope","Trump Proud")
xtable(therm.cor,digits=3)
dim(thermometer)

### SYMBOLIC IDEOLOGY AS A PREDICTOR IN 2016 ###
rm(list=ls())
additional<-read.table("anes16additional.txt",header=T)

#define model with default burn-in of 1000 iterations 
mod.ideol<-jags.model(file="MODELideology2016.txt", data=list(N=dim(additional)[1],
	vote.rep=additional$vote.rep,pid_x=additional$pid_x,issues=additional$issues,personal=additional$personal,
	ang.own=additional$ang.own,hp.own=additional$hp.own,afr.own=additional$afr.own,prd.own=additional$prd.own,
	disg.own=additional$disg.own,ideol=additional$ideol),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.ideol,n.iter=49000) 
#update(mod.ideol,n.iter=1000) 

#post burn-in sampling
output.mod.ideol<-coda.samples(mod.ideol,c("beta","constant","alpha"),500000)
#output.mod.ideol<-coda.samples(mod.ideol,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.ideol,10000)

#obtain model output
summary.cumulative.ideol<-summary(output.mod.ideol,quantiles=c(0.05,0.5,0.95))
summary.cumulative.ideol$statistics
summary.cumulative.ideol$quantiles
xtable(cbind(summary.cumulative.ideol$statistics[,1:2],summary.cumulative.ideol$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.ideol)
gelman.diag(output.mod.ideol)
Sys.time()



### DIRECTIONAL THEORY IN 2016 ###
rm(list=ls())
additional<-read.table("anes16additional.txt",header=T)

#define model with default burn-in of 1000 iterations 
mod.dir<-jags.model(file="MODELonly2016.txt", data=list(N=dim(additional)[1],
	vote.rep=additional$vote.rep,pid_x=additional$pid_x,issues=additional$issues.direction,personal=additional$personal,
	ang.own=additional$ang.own,hp.own=additional$hp.own,afr.own=additional$afr.own,prd.own=additional$prd.own,disg.own=additional$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.dir,n.iter=49000) 
#update(mod.dir,n.iter=1000) 

#post burn-in sampling
output.mod.dir<-coda.samples(mod.dir,c("beta","constant","alpha"),500000)
#output.mod.dir<-coda.samples(mod.dir,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.dir,10000)

#obtain model output
summary.cumulative.dir<-summary(output.mod.dir,quantiles=c(0.05,0.5,0.95))
summary.cumulative.dir$statistics
summary.cumulative.dir$quantiles
xtable(cbind(summary.cumulative.dir$statistics[,1:2],summary.cumulative.dir$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.dir)
gelman.diag(output.mod.dir)
Sys.time()



### QUADRATIC UTILITY IN 2016 ###
rm(list=ls())
additional<-read.table("anes16additional.txt",header=T)

#define model with default burn-in of 1000 iterations 
mod.sq<-jags.model(file="MODELonly2016.txt", data=list(N=dim(additional)[1],
	vote.rep=additional$vote.rep,pid_x=additional$pid_x,issues=additional$issues.sq,personal=additional$personal,
	ang.own=additional$ang.own,hp.own=additional$hp.own,afr.own=additional$afr.own,prd.own=additional$prd.own,disg.own=additional$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.sq,n.iter=49000) 
#update(mod.sq,n.iter=1000) 

#post burn-in sampling
output.mod.sq<-coda.samples(mod.sq,c("beta","constant","alpha"),500000)
#output.mod.sq<-coda.samples(mod.sq,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.sq,10000)

#obtain model output
summary.cumulative.sq<-summary(output.mod.sq,quantiles=c(0.05,0.5,0.95))
summary.cumulative.sq$statistics
summary.cumulative.sq$quantiles
xtable(cbind(summary.cumulative.sq$statistics[,1:2],summary.cumulative.sq$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.sq)
gelman.diag(output.mod.sq)
Sys.time()



### ABSOLUTE UTILITY AGAINST MEAN CANDIDATE PLACEMENT IN 2016 ###
rm(list=ls())
additional<-read.table("anes16additional.txt",header=T)

#define model with default burn-in of 1000 iterations 
mod.mean<-jags.model(file="MODELonly2016.txt", data=list(N=dim(additional)[1],
	vote.rep=additional$vote.rep,pid_x=additional$pid_x,issues=additional$issues.mean,personal=additional$personal,
	ang.own=additional$ang.own,hp.own=additional$hp.own,afr.own=additional$afr.own,prd.own=additional$prd.own,disg.own=additional$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.mean,n.iter=49000) 
#update(mod.mean,n.iter=1000) 

#post burn-in sampling
output.mod.mean<-coda.samples(mod.mean,c("beta","constant","alpha"),500000)
#output.mod.mean<-coda.samples(mod.mean,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.mean,10000)

#obtain model output
summary.cumulative.mean<-summary(output.mod.mean,quantiles=c(0.05,0.5,0.95))
summary.cumulative.mean$statistics
summary.cumulative.mean$quantiles
xtable(cbind(summary.cumulative.mean$statistics[,1:2],summary.cumulative.mean$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.mean)
gelman.diag(output.mod.mean)
Sys.time()



### PARTISAN EXTREMITY SUBGROUP MODELS IN 2016 ###
rm(list=ls())
data<-read.table("anes16.txt",header=T)

#creating subgroups
data.lean<-subset(data,subset=pid_x==0.333333333333333 | pid_x==0.666666666666667);dim(data.lean)
data.weak<-subset(data,subset=pid_x==0.166666666666667 | pid_x==0.833333333333333);dim(data.weak)
data.strong<-subset(data,subset=pid_x==0 | pid_x==1);dim(data.strong)

#recode party
data.lean$pid_x<-as.numeric(data.lean$pid_x>.5)
data.weak$pid_x<-as.numeric(data.weak$pid_x>.5)
data.strong$pid_x<-as.numeric(data.strong$pid_x>.5)

#Lean Model#
mod.lean<-jags.model(file="MODELonly2016.txt", data=list(N=dim(data.lean)[1],
	vote.rep=data.lean$vote.rep,pid_x=data.lean$pid_x,issues=data.lean$issues,personal=data.lean$personal,
	ang.own=data.lean$ang.own,hp.own=data.lean$hp.own,afr.own=data.lean$afr.own,prd.own=data.lean$prd.own,disg.own=data.lean$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.lean,n.iter=9000) 

#post burn-in sampling
output.mod.lean<-coda.samples(mod.lean,c("beta","constant","alpha"),50000)

#compute DIC
dic.samples(mod.lean,10000)

#obtain model output
summary.cumulative.lean<-summary(output.mod.lean,quantiles=c(0.05,0.5,0.95))
summary.cumulative.lean$statistics
summary.cumulative.lean$quantiles
xtable(cbind(summary.cumulative.lean$statistics[,1:2],summary.cumulative.lean$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.lean)
gelman.diag(output.mod.lean)
Sys.time()

#Weak Model#
mod.weak<-jags.model(file="MODELonly2016.txt", data=list(N=dim(data.weak)[1],
	vote.rep=data.weak$vote.rep,pid_x=data.weak$pid_x,issues=data.weak$issues,personal=data.weak$personal,
	ang.own=data.weak$ang.own,hp.own=data.weak$hp.own,afr.own=data.weak$afr.own,prd.own=data.weak$prd.own,disg.own=data.weak$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.weak,n.iter=9000) 

#post burn-in sampling
output.mod.weak<-coda.samples(mod.weak,c("beta","constant","alpha"),50000)

#compute DIC
dic.samples(mod.weak,10000)

#obtain model output
summary.cumulative.weak<-summary(output.mod.weak,quantiles=c(0.05,0.5,0.95))
summary.cumulative.weak$statistics
summary.cumulative.weak$quantiles
xtable(cbind(summary.cumulative.weak$statistics[,1:2],summary.cumulative.weak$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.weak)
gelman.diag(output.mod.weak)
Sys.time()

#Strong Model#
mod.strong<-jags.model(file="MODELonly2016.txt", data=list(N=dim(data.strong)[1],
	vote.rep=data.strong$vote.rep,pid_x=data.strong$pid_x,issues=data.strong$issues,personal=data.strong$personal,
	ang.own=data.strong$ang.own,hp.own=data.strong$hp.own,afr.own=data.strong$afr.own,prd.own=data.strong$prd.own,disg.own=data.strong$disg.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.strong,n.iter=9000) 

#post burn-in sampling
output.mod.strong<-coda.samples(mod.strong,c("beta","constant","alpha"),50000)

#compute DIC
dic.samples(mod.strong,10000)

#obtain model output
summary.cumulative.strong<-summary(output.mod.strong,quantiles=c(0.05,0.5,0.95))
summary.cumulative.strong$statistics
summary.cumulative.strong$quantiles
xtable(cbind(summary.cumulative.strong$statistics[,1:2],summary.cumulative.strong$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.strong)
gelman.diag(output.mod.strong)
Sys.time()



### 2012 MODEL ###
rm(list=ls())
data.2012<-read.table("anes12.txt",header=T)

#descriptives
data.desc<-subset(data.2012,select=-caseid)
print(cbind(
apply(data.desc,2,mean),
apply(data.desc,2,sd),
apply(data.desc,2,min),
apply(data.desc,2,max)
),digits=3)

#define model with default burn-in of 1000 iterations 
mod.12<-jags.model(file="MODELannual.txt", data=list(N=dim(data.2012)[1],
	vote.rep=data.2012$vote.rep,pid_x=data.2012$pid_x,issues=data.2012$issues,personal=data.2012$personal,
	ang.own=data.2012$ang.own,hp.own=data.2012$hp.own,afr.own=data.2012$afr.own,prd.own=data.2012$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.12,n.iter=49000) 
#update(mod.12,n.iter=1000) 

#post burn-in sampling
output.mod.12<-coda.samples(mod.12,c("beta","constant","alpha"),500000)
#output.mod.12<-coda.samples(mod.12,c("beta","constant","alpha"),1000)

#Density Plots
pdf("anxPid2012.pdf",pointsize=8,width=4,height=4)
anx.pid<-c(output.mod.12[[1]][,"beta[5]"],output.mod.12[[2]][,"beta[5]"],output.mod.12[[3]][,"beta[5]"])
plot(density(anx.pid),xlab="Coefficient",main="Interaction Between Anxiety and Partisanship")
dev.off()

pdf("anxIss2012.pdf",pointsize=8,width=4,height=4)
anx.iss<-c(output.mod.12[[1]][,"beta[6]"],output.mod.12[[2]][,"beta[6]"],output.mod.12[[3]][,"beta[6]"])
plot(density(anx.iss),xlab="Coefficient",main="Interaction Between Anxiety and Issue Distance")
dev.off()

pdf("anxPers2012.pdf",pointsize=8,width=4,height=4)
anx.pers<-c(output.mod.12[[1]][,"beta[7]"],output.mod.12[[2]][,"beta[7]"],output.mod.12[[3]][,"beta[7]"])
plot(density(anx.pers),xlab="Coefficient",main="Interaction Between Anxiety and Candidate Qualities")
dev.off()

#compute DIC
dic.samples(mod.12,10000)

#obtain model output
summary.cumulative<-summary(output.mod.12,quantiles=c(0.05,0.5,0.95))
summary.cumulative$statistics
summary.cumulative$quantiles
xtable(cbind(summary.cumulative$statistics[,1:2],summary.cumulative$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.12)
gelman.diag(output.mod.12)
Sys.time()

###Predicted Probabilities 2012###
#Anxiety is on a logit 0 to 1 scale
plogis(0)
plogis(1)
plogis(-1)

#bind samples into a matrix
output.mat<-rbind(output.mod.12[[1]],output.mod.12[[2]],output.mod.12[[3]])
output.mat<-output.mat[,-c(1:4)]
output.means<-apply(output.mat,2,mean)

#create vectors of average predictor values at low, medium, and high anxiety
low.averages<-c(0.4391,0.5009,0.6259,0.2689414,0.2689414*0.4391,0.2689414*0.5009,0.2689414*0.6259,1)
med.averages<-c(0.4391,0.5009,0.6259,0.5,0.5*0.4391,0.5*0.5009,0.5*0.6259,1)
high.averages<-c(0.4391,0.5009,0.6259,0.7310586,0.7310586*0.4391,0.7310586*0.5009,0.7310586*0.6259,1)
###plogis(low.averages%*%output.means)
###plogis(med.averages%*%output.means)
###plogis(high.averages%*%output.means)

##create matrices where PARTY varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.party<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.party[,1]<-ruler
low.averages.party[,5]<-ruler*0.2689414
med.averages.party<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.party[,1]<-ruler
med.averages.party[,5]<-ruler*0.5
high.averages.party<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.party[,1]<-ruler
high.averages.party[,5]<-ruler*0.7310586

#get averages
low.averages.party.prob<-plogis(low.averages.party%*%output.means)
med.averages.party.prob<-plogis(med.averages.party%*%output.means)
high.averages.party.prob<-plogis(high.averages.party%*%output.means)

#forecast for each iteration
low.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.party.all[i,]<-plogis(low.averages.party%*%output.mat[i,])
	}
low.party.ci.low<-apply(low.party.all,2,quantile,.05)
low.party.ci.high<-apply(low.party.all,2,quantile,.95)

med.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.party.all[i,]<-plogis(med.averages.party%*%output.mat[i,])
	}
med.party.ci.low<-apply(med.party.all,2,quantile,.05)
med.party.ci.high<-apply(med.party.all,2,quantile,.95)

high.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.party.all[i,]<-plogis(high.averages.party%*%output.mat[i,])
	}
high.party.ci.low<-apply(high.party.all,2,quantile,.05)
high.party.ci.high<-apply(high.party.all,2,quantile,.95)

pdf("partyProb2012.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.party.ci.low,rev(med.party.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.party.ci.low,rev(low.party.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.party.ci.low,rev(high.party.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("partyProb2012mean.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where ISSUE proximity varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.issue<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.issue[,2]<-ruler
low.averages.issue[,6]<-ruler*0.2689414
med.averages.issue<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.issue[,2]<-ruler
med.averages.issue[,6]<-ruler*0.5
high.averages.issue<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.issue[,2]<-ruler
high.averages.issue[,6]<-ruler*0.7310586

#get averages
low.averages.issue.prob<-plogis(low.averages.issue%*%output.means)
med.averages.issue.prob<-plogis(med.averages.issue%*%output.means)
high.averages.issue.prob<-plogis(high.averages.issue%*%output.means)

#forecast for each iteration
low.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.issue.all[i,]<-plogis(low.averages.issue%*%output.mat[i,])
	}
low.issue.ci.low<-apply(low.issue.all,2,quantile,.05)
low.issue.ci.high<-apply(low.issue.all,2,quantile,.95)

med.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.issue.all[i,]<-plogis(med.averages.issue%*%output.mat[i,])
	}
med.issue.ci.low<-apply(med.issue.all,2,quantile,.05)
med.issue.ci.high<-apply(med.issue.all,2,quantile,.95)

high.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.issue.all[i,]<-plogis(high.averages.issue%*%output.mat[i,])
	}
high.issue.ci.low<-apply(high.issue.all,2,quantile,.05)
high.issue.ci.high<-apply(high.issue.all,2,quantile,.95)

pdf("issueProb2012.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Romney Issue Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.issue.ci.low,rev(med.issue.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.issue.ci.low,rev(low.issue.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.issue.ci.low,rev(high.issue.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("issueProb2012mean.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Romney Issue Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where PERSONAL qualities vary for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.personal<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.personal[,3]<-ruler
low.averages.personal[,7]<-ruler*0.2689414
med.averages.personal<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.personal[,3]<-ruler
med.averages.personal[,7]<-ruler*0.5
high.averages.personal<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.personal[,3]<-ruler
high.averages.personal[,7]<-ruler*0.7310586

#get averages
low.averages.personal.prob<-plogis(low.averages.personal%*%output.means)
med.averages.personal.prob<-plogis(med.averages.personal%*%output.means)
high.averages.personal.prob<-plogis(high.averages.personal%*%output.means)

#forecast for each iteration
low.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.personal.all[i,]<-plogis(low.averages.personal%*%output.mat[i,])
	}
low.personal.ci.low<-apply(low.personal.all,2,quantile,.05)
low.personal.ci.high<-apply(low.personal.all,2,quantile,.95)

med.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.personal.all[i,]<-plogis(med.averages.personal%*%output.mat[i,])
	}
med.personal.ci.low<-apply(med.personal.all,2,quantile,.05)
med.personal.ci.high<-apply(med.personal.all,2,quantile,.95)

high.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.personal.all[i,]<-plogis(high.averages.personal%*%output.mat[i,])
	}
high.personal.ci.low<-apply(high.personal.all,2,quantile,.05)
high.personal.ci.high<-apply(high.personal.all,2,quantile,.95)

pdf("personalProb2012.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Romney Personal Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.personal.ci.low,rev(med.personal.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.personal.ci.low,rev(low.personal.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.personal.ci.low,rev(high.personal.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("personalProb2012mean.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Romney Personal Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()



###Alternate Predicted Probabilities 2012, computed at EMPIRICAL minimum, mean, and maximum.###
#Anxiety is on a logit 0 to 1 scale
plogis(0)#now 0.4905
plogis(1)#now 0.9128
plogis(-1)#now 0.1977

#bind samples into a matrix
output.mat<-rbind(output.mod.12[[1]],output.mod.12[[2]],output.mod.12[[3]])
output.mat<-output.mat[,-c(1:4)]
output.means<-apply(output.mat,2,mean)

#create vectors of average predictor values at low, medium, and high anxiety
low.averages<-c(0.4391,0.5009,0.6259,0.1977,0.1977*0.4391,0.1977*0.5009,0.1977*0.6259,1)
med.averages<-c(0.4391,0.5009,0.6259, 0.4905, 0.4905*0.4391, 0.4905*0.5009, 0.4905*0.6259,1)
high.averages<-c(0.4391,0.5009,0.6259,0.9128,0.9128*0.4391,0.9128*0.5009,0.9128*0.6259,1)
###plogis(low.averages%*%output.means)
###plogis(med.averages%*%output.means)
###plogis(high.averages%*%output.means)

##create matrices where PARTY varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.party<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.party[,1]<-ruler
low.averages.party[,5]<-ruler*0.1977
med.averages.party<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.party[,1]<-ruler
med.averages.party[,5]<-ruler* 0.4905
high.averages.party<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.party[,1]<-ruler
high.averages.party[,5]<-ruler*0.9128

#get averages
low.averages.party.prob<-plogis(low.averages.party%*%output.means)
med.averages.party.prob<-plogis(med.averages.party%*%output.means)
high.averages.party.prob<-plogis(high.averages.party%*%output.means)

#forecast for each iteration
low.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.party.all[i,]<-plogis(low.averages.party%*%output.mat[i,])
	}
low.party.ci.low<-apply(low.party.all,2,quantile,.05)
low.party.ci.high<-apply(low.party.all,2,quantile,.95)

med.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.party.all[i,]<-plogis(med.averages.party%*%output.mat[i,])
	}
med.party.ci.low<-apply(med.party.all,2,quantile,.05)
med.party.ci.high<-apply(med.party.all,2,quantile,.95)

high.party.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.party.all[i,]<-plogis(high.averages.party%*%output.mat[i,])
	}
high.party.ci.low<-apply(high.party.all,2,quantile,.05)
high.party.ci.high<-apply(high.party.all,2,quantile,.95)

pdf("partyProb2012Alt.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.party.ci.low,rev(med.party.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.party.ci.low,rev(low.party.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.party.ci.low,rev(high.party.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("partyProb2012meanAlt.pdf")
plot(y=med.averages.party.prob,x=ruler,type='l',xlab="Partisanship--Higher Values Republican",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.party.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.party.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where ISSUE proximity varies for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.issue<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.issue[,2]<-ruler
low.averages.issue[,6]<-ruler*0.1977
med.averages.issue<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.issue[,2]<-ruler
med.averages.issue[,6]<-ruler* 0.4905
high.averages.issue<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.issue[,2]<-ruler
high.averages.issue[,6]<-ruler*0.9128

#get averages
low.averages.issue.prob<-plogis(low.averages.issue%*%output.means)
med.averages.issue.prob<-plogis(med.averages.issue%*%output.means)
high.averages.issue.prob<-plogis(high.averages.issue%*%output.means)

#forecast for each iteration
low.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.issue.all[i,]<-plogis(low.averages.issue%*%output.mat[i,])
	}
low.issue.ci.low<-apply(low.issue.all,2,quantile,.05)
low.issue.ci.high<-apply(low.issue.all,2,quantile,.95)

med.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.issue.all[i,]<-plogis(med.averages.issue%*%output.mat[i,])
	}
med.issue.ci.low<-apply(med.issue.all,2,quantile,.05)
med.issue.ci.high<-apply(med.issue.all,2,quantile,.95)

high.issue.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.issue.all[i,]<-plogis(high.averages.issue%*%output.mat[i,])
	}
high.issue.ci.low<-apply(high.issue.all,2,quantile,.05)
high.issue.ci.high<-apply(high.issue.all,2,quantile,.95)

pdf("issueProb2012Alt.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Romney Issue Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.issue.ci.low,rev(med.issue.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.issue.ci.low,rev(low.issue.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.issue.ci.low,rev(high.issue.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("issueProb2012meanAlt.pdf")
plot(y=med.averages.issue.prob,x=ruler,type='l',xlab="Romney Issue Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.issue.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.issue.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

##create matrices where PERSONAL qualities vary for low, medium, and high anxiety
ruler<-seq(0,1,.01)
low.averages.personal<-matrix(rep(low.averages,length(ruler)),nrow=length(ruler),ncol=length(low.averages),byrow=T)
low.averages.personal[,3]<-ruler
low.averages.personal[,7]<-ruler*0.1977
med.averages.personal<-matrix(rep(med.averages,length(ruler)),nrow=length(ruler),ncol=length(med.averages),byrow=T)
med.averages.personal[,3]<-ruler
med.averages.personal[,7]<-ruler* 0.4905
high.averages.personal<-matrix(rep(high.averages,length(ruler)),nrow=length(ruler),ncol=length(high.averages),byrow=T)
high.averages.personal[,3]<-ruler
high.averages.personal[,7]<-ruler*0.9128

#get averages
low.averages.personal.prob<-plogis(low.averages.personal%*%output.means)
med.averages.personal.prob<-plogis(med.averages.personal%*%output.means)
high.averages.personal.prob<-plogis(high.averages.personal%*%output.means)

#forecast for each iteration
low.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	low.personal.all[i,]<-plogis(low.averages.personal%*%output.mat[i,])
	}
low.personal.ci.low<-apply(low.personal.all,2,quantile,.05)
low.personal.ci.high<-apply(low.personal.all,2,quantile,.95)

med.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	med.personal.all[i,]<-plogis(med.averages.personal%*%output.mat[i,])
	}
med.personal.ci.low<-apply(med.personal.all,2,quantile,.05)
med.personal.ci.high<-apply(med.personal.all,2,quantile,.95)

high.personal.all<-matrix(NA,nrow=dim(output.mat)[1],ncol=length(ruler))
for(i in 1:dim(output.mat)[1]){
	high.personal.all[i,]<-plogis(high.averages.personal%*%output.mat[i,])
	}
high.personal.ci.low<-apply(high.personal.all,2,quantile,.05)
high.personal.ci.high<-apply(high.personal.all,2,quantile,.95)

pdf("personalProb2012Alt.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Romney Personal Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(med.personal.ci.low,rev(med.personal.ci.high)),border=F,col=rgb(.6,.6,.6,.3))
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(low.personal.ci.low,rev(low.personal.ci.high)),border=F,col=rgb(.9,0,0,.3))
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
polygon(x=c(ruler,rev(ruler)),y=c(high.personal.ci.low,rev(high.personal.ci.high)),border=F,col=rgb(0,0,.9,.3))
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()

pdf("personalProb2012meanAlt.pdf")
plot(y=med.averages.personal.prob,x=ruler,type='l',xlab="Romney Personal Advantage",ylab="Probability of Romney Vote",ylim=c(0,1),lwd=2)
lines(y=low.averages.personal.prob,x=ruler,lty=2,col="red",lwd=2)
lines(y=high.averages.personal.prob,x=ruler,lty=3,col="blue",lwd=2)
legend(x=0,y=1,legend=c("Low Anxiety","Medium Anxiety","High Anxiety"),lwd=2,lty=c(2,1,3),col=c("red","black","blue"))
dev.off()











### POOLING 2012 & 2016 to get anxiety on the same scale ###
rm(list=ls())
data<-read.table("anes16.txt",header=T)
data.2012<-read.table("anes12.txt",header=T)

#combine data
data<-subset(data,select=-c(disg.own))
data$year<-2016
data.2012$year<-2012
pooled<-rbind(data,data.2012)

#define model with default burn-in of 1000 iterations 
mod.pool<-jags.model(file="MODELannual.txt", data=list(N=dim(pooled)[1],
	vote.rep=pooled$vote.rep,pid_x=pooled$pid_x,issues=pooled$issues,personal=pooled$personal,
	ang.own=pooled$ang.own,hp.own=pooled$hp.own,afr.own=pooled$afr.own,prd.own=pooled$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.pool,n.iter=49000) 
#update(mod.pool,n.iter=1000) 

#post burn-in sampling
output.mod.pool<-coda.samples(mod.pool,c("anxiety.own"),50000)
#output.mod.pool<-coda.samples(mod.pool,c("anxiety.own"),1000)

#obtain model output
summary.cumulative.pool<-summary(output.mod.pool,quantiles=c(0.5))
pooled$anxiety.own<-summary.cumulative.pool$statistics[,1]
#head(pooled)
#head(summary.cumulative.pool$statistics)
write.csv(pooled,"anxiety1216.csv",row.names=F)
Sys.time()



### ANNUAL MODELS ###
rm(list=ls())
data.cum<-read.table("anesCumulative.txt",header=T)
ang.mod<-CrossTable(x=data.cum$ang.own,y=data.cum$year,prop.r=F,prop.t=F,chisq=T)
afr.mod<-CrossTable(x=data.cum$afr.own,y=data.cum$year,prop.r=F,prop.t=F,chisq=T)
hp.mod<-CrossTable(x=data.cum$hp.own,y=data.cum$year,prop.r=F,prop.t=F,chisq=T)
prd.mod<-CrossTable(x=data.cum$prd.own,y=data.cum$year,prop.r=F,prop.t=F,chisq=T)
#table(as.numeric(data$ang.own>0))/sum(table(as.numeric(data$ang.own>0)))
#table(as.numeric(data.2012$ang.own>0))/sum(table(as.numeric(data.2012$ang.own>0)))

#pdf("emotionsTime.pdf",width=4,height=4,pointsize=8)
plot(y=ang.mod$prop.col[2,],x=names(ang.mod$prop.col[2,]),type='l',ylim=c(0,1),lwd=2,xlab="Year",ylab="Proportion Reporting Emotion",axes=F)
axis(1,at=seq(1984,2016,4));axis(2);box()
lines(y=afr.mod$prop.col[2,],x=names(afr.mod$prop.col[2,]),lwd=2,lty=2,col="red")
lines(y=hp.mod$prop.col[2,],x=names(hp.mod$prop.col[2,]),lwd=2,lty=3,col="blue")
lines(y=prd.mod$prop.col[2,],x=names(prd.mod$prop.col[2,]),lwd=2,lty=4,col="ForestGreen")
text(2008,0,"Afraid")
arrows(x0=2008,y0=0.01,x1=2009,y1=.1,length=.1)
text(2008,.7,"Proud")
arrows(x0=2008,y0=0.71,x1=2009,y1=.8,length=.1)
text(1985,0.4,"Angry")
arrows(x0=1985,y0=.38,x1=1986,y1=.23,length=.1)
text(1985,1,"Hope")
arrows(x0=1985,y0=.98,x1=1986,y1=.8,length=.1)
#dev.off()

#data subsets
year.1984<-subset(data.cum,subset=year==1984);dim(year.1984)
year.1988<-subset(data.cum,subset=year==1988);dim(year.1988)
year.1992<-subset(data.cum,subset=year==1992);dim(year.1992)
year.1996<-subset(data.cum,subset=year==1996);dim(year.1996)
year.2000<-subset(data.cum,subset=year==2000);dim(year.2000)
year.2004<-subset(data.cum,subset=year==2004);dim(year.2004)
year.2008<-subset(data.cum,subset=year==2008);dim(year.2008)
year.2012<-subset(data.cum,subset=year==2012);dim(year.2012)
year.2016<-subset(data.cum,subset=year==2016);dim(year.2016)

#1984 model
#define model with default burn-in of 1000 iterations 
mod.y.84<-jags.model(file="MODELannual.txt", data=list(N=dim(year.1984)[1],
	vote.rep=year.1984$vote.rep,pid_x=year.1984$pid_x,issues=year.1984$issues,personal=year.1984$personal,
	ang.own=year.1984$ang.own,hp.own=year.1984$hp.own,afr.own=year.1984$afr.own,prd.own=year.1984$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.84,n.iter=9000) 
#update(mod.y.84,n.iter=1000) 

#post burn-in sampling
output.mod.y.84<-coda.samples(mod.y.84,c("beta","constant","alpha"),50000)
#output.mod.y.84<-coda.samples(mod.y.84,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.84,10000)

#obtain model output
summary.cumulative.84<-summary(output.mod.y.84,quantiles=c(0.05,0.5,0.95))
summary.cumulative.84$statistics
summary.cumulative.84$quantiles
xtable(cbind(summary.cumulative.84$statistics[,1:2],summary.cumulative.84$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.84)
gelman.diag(output.mod.y.84)
Sys.time()

#1988 model
#define model with default burn-in of 1000 iterations 
mod.y.88<-jags.model(file="MODELannual.txt", data=list(N=dim(year.1988)[1],
	vote.rep=year.1988$vote.rep,pid_x=year.1988$pid_x,issues=year.1988$issues,personal=year.1988$personal,
	ang.own=year.1988$ang.own,hp.own=year.1988$hp.own,afr.own=year.1988$afr.own,prd.own=year.1988$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.88,n.iter=9000) 
#update(mod.y.88,n.iter=1000) 

#post burn-in sampling
output.mod.y.88<-coda.samples(mod.y.88,c("beta","constant","alpha"),50000)
#output.mod.y.88<-coda.samples(mod.y.88,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.88,10000)

#obtain model output
summary.cumulative.88<-summary(output.mod.y.88,quantiles=c(0.05,0.5,0.95))
summary.cumulative.88$statistics
summary.cumulative.88$quantiles
xtable(cbind(summary.cumulative.88$statistics[,1:2],summary.cumulative.88$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.88)
gelman.diag(output.mod.y.88)
Sys.time()

#1992 model
#define model with default burn-in of 1000 iterations 
mod.y.92<-jags.model(file="MODELannual.txt", data=list(N=dim(year.1992)[1],
	vote.rep=year.1992$vote.rep,pid_x=year.1992$pid_x,issues=year.1992$issues,personal=year.1992$personal,
	ang.own=year.1992$ang.own,hp.own=year.1992$hp.own,afr.own=year.1992$afr.own,prd.own=year.1992$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.92,n.iter=9000) 
#update(mod.y.92,n.iter=1000) 

#post burn-in sampling
output.mod.y.92<-coda.samples(mod.y.92,c("beta","constant","alpha"),50000)
#output.mod.y.92<-coda.samples(mod.y.92,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.92,10000)

#obtain model output
summary.cumulative.92<-summary(output.mod.y.92,quantiles=c(0.05,0.5,0.95))
summary.cumulative.92$statistics
summary.cumulative.92$quantiles
xtable(cbind(summary.cumulative.92$statistics[,1:2],summary.cumulative.92$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.92)
gelman.diag(output.mod.y.92)
Sys.time()

#1996 model
#define model with default burn-in of 1000 iterations 
mod.y.96<-jags.model(file="MODELannual.txt", data=list(N=dim(year.1996)[1],
	vote.rep=year.1996$vote.rep,pid_x=year.1996$pid_x,issues=year.1996$issues,personal=year.1996$personal,
	ang.own=year.1996$ang.own,hp.own=year.1996$hp.own,afr.own=year.1996$afr.own,prd.own=year.1996$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.96,n.iter=9000) 
#update(mod.y.96,n.iter=1000) 

#post burn-in sampling
output.mod.y.96<-coda.samples(mod.y.96,c("beta","constant","alpha"),50000)
#output.mod.y.96<-coda.samples(mod.y.96,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.96,10000)

#obtain model output
summary.cumulative.96<-summary(output.mod.y.96,quantiles=c(0.05,0.5,0.95))
summary.cumulative.96$statistics
summary.cumulative.96$quantiles
xtable(cbind(summary.cumulative.96$statistics[,1:2],summary.cumulative.96$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.96)
gelman.diag(output.mod.y.96)
Sys.time()

#2000 model
#define model with default burn-in of 1000 iterations 
mod.y.00<-jags.model(file="MODELannual.txt", data=list(N=dim(year.2000)[1],
	vote.rep=year.2000$vote.rep,pid_x=year.2000$pid_x,issues=year.2000$issues,personal=year.2000$personal,
	ang.own=year.2000$ang.own,hp.own=year.2000$hp.own,afr.own=year.2000$afr.own,prd.own=year.2000$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.00,n.iter=9000) 
#update(mod.y.00,n.iter=1000) 

#post burn-in sampling
output.mod.y.00<-coda.samples(mod.y.00,c("beta","constant","alpha"),50000)
#output.mod.y.00<-coda.samples(mod.y.00,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.00,10000)

#obtain model output
summary.cumulative.00<-summary(output.mod.y.00,quantiles=c(0.05,0.5,0.95))
summary.cumulative.00$statistics
summary.cumulative.00$quantiles
xtable(cbind(summary.cumulative.00$statistics[,1:2],summary.cumulative.00$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.00)
gelman.diag(output.mod.y.00)
Sys.time()

#2004 model
#define model with default burn-in of 1000 iterations 
mod.y.04<-jags.model(file="MODELannual.txt", data=list(N=dim(year.2004)[1],
	vote.rep=year.2004$vote.rep,pid_x=year.2004$pid_x,issues=year.2004$issues,personal=year.2004$personal,
	ang.own=year.2004$ang.own,hp.own=year.2004$hp.own,afr.own=year.2004$afr.own,prd.own=year.2004$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.04,n.iter=9000) 
#update(mod.y.04,n.iter=1000) 

#post burn-in sampling
output.mod.y.04<-coda.samples(mod.y.04,c("beta","constant","alpha"),50000)
#output.mod.y.04<-coda.samples(mod.y.04,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.04,10000)

#obtain model output
summary.cumulative.04<-summary(output.mod.y.04,quantiles=c(0.05,0.5,0.95))
summary.cumulative.04$statistics
summary.cumulative.04$quantiles
xtable(cbind(summary.cumulative.04$statistics[,1:2],summary.cumulative.04$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.04)
gelman.diag(output.mod.y.04)
Sys.time()

#2008 model
#define model with default burn-in of 1000 iterations 
mod.y.08<-jags.model(file="MODELannual.txt", data=list(N=dim(year.2008)[1],
	vote.rep=year.2008$vote.rep,pid_x=year.2008$pid_x,issues=year.2008$issues,personal=year.2008$personal,
	ang.own=year.2008$ang.own,hp.own=year.2008$hp.own,afr.own=year.2008$afr.own,prd.own=year.2008$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.08,n.iter=9000) 
#update(mod.y.08,n.iter=1000) 

#post burn-in sampling
output.mod.y.08<-coda.samples(mod.y.08,c("beta","constant","alpha"),50000)
#output.mod.y.08<-coda.samples(mod.y.08,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.08,10000)

#obtain model output
summary.cumulative.08<-summary(output.mod.y.08,quantiles=c(0.05,0.5,0.95))
summary.cumulative.08$statistics
summary.cumulative.08$quantiles
xtable(cbind(summary.cumulative.08$statistics[,1:2],summary.cumulative.08$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.08)
gelman.diag(output.mod.y.08)
Sys.time()

#2012 model
#define model with default burn-in of 1000 iterations 
mod.y.12<-jags.model(file="MODELannual.txt", data=list(N=dim(year.2012)[1],
	vote.rep=year.2012$vote.rep,pid_x=year.2012$pid_x,issues=year.2012$issues,personal=year.2012$personal,
	ang.own=year.2012$ang.own,hp.own=year.2012$hp.own,afr.own=year.2012$afr.own,prd.own=year.2012$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.12,n.iter=9000) 
#update(mod.y.12,n.iter=1000) 

#post burn-in sampling
output.mod.y.12<-coda.samples(mod.y.12,c("beta","constant","alpha"),50000)
#output.mod.y.12<-coda.samples(mod.y.12,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.12,10000)

#obtain model output
summary.cumulative.12<-summary(output.mod.y.12,quantiles=c(0.05,0.5,0.95))
summary.cumulative.12$statistics
summary.cumulative.12$quantiles
xtable(cbind(summary.cumulative.12$statistics[,1:2],summary.cumulative.12$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.12)
gelman.diag(output.mod.y.12)
Sys.time()

#2016 model
#define model with default burn-in of 1000 iterations 
mod.y.16<-jags.model(file="MODELannual.txt", data=list(N=dim(year.2016)[1],
	vote.rep=year.2016$vote.rep,pid_x=year.2016$pid_x,issues=year.2016$issues,personal=year.2016$personal,
	ang.own=year.2016$ang.own,hp.own=year.2016$hp.own,afr.own=year.2016$afr.own,prd.own=year.2016$prd.own),
	n.chains=3, n.adapt=1000) #initial burn-in of 1000

#burn in
update(mod.y.16,n.iter=9000) 
#update(mod.y.16,n.iter=1000) 

#post burn-in sampling
output.mod.y.16<-coda.samples(mod.y.16,c("beta","constant","alpha"),50000)
#output.mod.y.16<-coda.samples(mod.y.16,c("beta","constant","alpha"),1000)

#compute DIC
dic.samples(mod.y.16,10000)

#obtain model output
summary.cumulative.16<-summary(output.mod.y.16,quantiles=c(0.05,0.5,0.95))
summary.cumulative.16$statistics
summary.cumulative.16$quantiles
xtable(cbind(summary.cumulative.16$statistics[,1:2],summary.cumulative.16$quantiles[,c(1,3)]),digits=3)
geweke.diag(output.mod.y.16)
gelman.diag(output.mod.y.16)
Sys.time()


### Cronbach's alpha for 2016 ###
library(psych)
data<-read.table("anes16.txt",header=T)
emotions<-subset(data,select=c(ang.own,hp.own,afr.own,prd.own,disg.own))
emotions$hp.own<- -1*emotions$hp.own
emotions$prd.own<- -1*emotions$prd.own
alpha(emotions)


### Joint distribution of anxiety and issue preference ###
joint.anx<-read.csv("anxiety1216.csv",header=T)
joint.anx.12<-subset(joint.anx,year==2012)
joint.anx.16<-subset(joint.anx,year==2016)

#pdf("anxIssuePool.pdf")
cor(joint.anx$anxiety.own,joint.anx$issues)
plot(x=joint.anx$anxiety.own,y=joint.anx$issues,ylim=c(0,1),xlim=c(.2,.9),xlab="Anxiety",ylab="Issue Proximity",main="Pooled")
abline(lm(issues~anxiety.own,data=joint.anx),lwd=2)
#dev.off()

#pdf("anxIssue2012.pdf")
cor(joint.anx.12$anxiety.own,joint.anx.12$issues)
plot(x=joint.anx.12$anxiety.own,y=joint.anx.12$issues,ylim=c(0,1),xlim=c(.2,.9),xlab="Anxiety",ylab="Issue Proximity",main="2012")
abline(lm(issues~anxiety.own,data=joint.anx.12),lwd=2)
#dev.off()

#pdf("anxIssue2016.pdf")
cor(joint.anx.16$anxiety.own,joint.anx.16$issues)
plot(x=joint.anx.16$anxiety.own,y=joint.anx.16$issues,ylim=c(0,1),xlim=c(.2,.9),xlab="Anxiety",ylab="Issue Proximity",main="2016")
abline(lm(issues~anxiety.own,data=joint.anx.16),lwd=2)
#dev.off()
